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Abstract 

We describe a generalized scheme for the probability-changing cluster (PCC) algo- 
rithm, based on the study of the finite-size scaling property of the correlation ratio, 
the ratio of the correlation functions with different distances. We apply this gener- 
alized PCC algorithm to the two-dimensional 6-state clock model. We also discuss 
the combination of the cluster algorithm and the extended ensemble method. We 
derive a rigorous broad histogram relation for the bond number. A Monte Carlo dy- 
namics based on the number of potential moves for the bond number is proposed, 
and applied to the three-dimensional Ising and 3-state Potts models. 

Key words: Cluster algorithm; Finite-size scaling; Correlation ratio; Broad 
histogram relation 



1 Introduction 



In the Monte Carlo simulation, we sometimes suffer from the problem of slow 
dynamics. The critical slowing down near the critical point, the phase separa- 
tion dynamics at low temperature, the slow dynamics due to the randomness 
or frustration, and the low-temperature slow dynamics in quantum Monte 
Carlo simulation are examples of the problems of slow dynamics. 

We may classify the attempts to conquer the slow dynamics in the Monte Carlo 
simulation into two categories. The first category is the cluster algorithm, such 
as the Swendsen-Wang (SW) algorithm [1] and the Wolff algorithm [2]. The 
second one is the extended ensemble method. The multicanonical method 
[3,4], the broad histogram method [5], and the flat histogram method [6,7] are 
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examples of the second category. Recently, Wang and Landau [8] proposed an 
efficient algorithm to accelerate the calculation of the energy density of states 
(DOS). Yamaguchi and Okabe [9] have successfully used the Wang-Landau 
algorithm for the study of the antiferromagnetic g-state Potts models. 

Tomita and Okabe [10] recently proposed an effective cluster algorithm, which 
is called the probability-changing cluster (PCC) algorithm, of tuning the crit- 
ical point automatically. The PCC algorithm is an extension of the SW al- 
gorithm [1]; we change the probability of connecting spins of the same type, 
essentially the temperature, in the process of the Monte Carlo spin update. 
We showed the effectiveness of the PCC algorithm for the two-dimensional 
(2D) and three-dimensional (3D) Potts models [10], determining the critical 
point and exponents. We can extract information on critical phenomena with 
much less numerical effort. The PCC algorithm is quite useful for studying 
random systems, where the distribution of the critical temperature, T c , is im- 
portant. We applied the PCC algorithm to the 2D diluted Ising model [11], 
investigating the crossover and self-averaging properties of random systems. 
We also extended the PCC algorithm to the problem of the vector order pa- 
rameter [12]; studying the 2D XY and clock models, we showed that the PCC 
algorithm is also useful for the Kosterlitz-Thouless (KT) transition [13]. 

The combination of approaches of two categories, the cluster algorithm and 
the extended ensemble method, is a challenging problem. Janke and Kappler 
[14] proposed a trial to combine the multicanonical method and the cluster 
algorithm; their method is called the multibondic ensemble method. Recently, 
Yamaguchi and Kawashima [15] improved the multibondic ensemble method, 
and also showed that the combination of the Wang-Landau algorithm and the 
improved multibondic ensemble method yields much better statistics com- 
pared to the original multibondic ensemble method [14]. 

Here, we pick up two recent topics of new Monte Carlo algorithms. We first 
discuss the generalization of the PCC algorithm. This generalized scheme is 
based on the finite-size scaling (FSS) property of the correlation ratio. Second, 
for the algorithm to combine the cluster algorithm and the extended ensemble 
method, we derive a rigorous broad histogram relation for the bond number, 
and propose the flat histogram method for the bond number. 



2 Generalized scheme for the PCC algorithm 

We start with reviewing the idea of the PCC algorithm [10]. We explain the 
case of the ferromagnetic g-state Potts model, whose Hamiltonian is given by 

-J ^2(5^-1), a t = l,2,---,q. (1) 

<i,j> 
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We construct the Kasteleyn and Fortuin (KF) [16] clusters using the prob- 
ability p, where p is the probability of connecting spins of the same type, 
p — 1 — e~ J / k B T , The correspondence of the spontaneous magnetization of the 
g-state Potts model and the percolation probability of the bond percolation 
model was discussed by Hu [17]. Then, we check whether the system is per- 
colating or not. If the system is percolating (not percolating) in the previous 
test, we decrease (increase) p by Ap (> 0). Spins are updated following the 
same rule as the SW algorithm. After repeating the above processes, the dis- 
tribution of p for Monte Carlo samples approaches the Gaussian distribution 
of which mean value is p c (L); p c (L) is the probability of connecting spins, such 
that the existence probability E p becomes 1/2. The existence probability E p 
[18] is the probability that the system percolates. Since E p follows the FSS 
near the critical point, 

E p (p,L)^X(tL^), t = (p c -p)/p c , (2) 



where p c is the critical value of p for the infinite system and v is the correlation- 
length critical exponent, we can estimate p c from the size dependence of p c {L) 
using Eq. (2) and, in turn, estimate T c through the relation p c = 1 — e ~ J / fcsTc . 

In the original formulation of the PCC algorithm [10], we used the KF repre- 
sentation in two ways. First, we make a cluster flip as in the SW algorithm 
[1]. Second, we change the probability of connecting spins of the same type, 
p, depending on the observation whether clusters are percolating or not. The 
point is that E p has the FSS property with a single scaling variable. We may 
use quantities other than E v which have a similar FSS relation. In the FSS 
analysis of the Monte Carlo simulation, we often use the Binder ratio [19], 
which is essentially the ratio of the moments of the order parameter m. The 
moment ratio has the FSS property with a single scaling variable, 



as far as the corrections to FSS are negligible. The angular brackets indicate 
the thermal average. The moment ratio derived from a snapshot spin configu- 
ration is always one; therefore, the instantaneous moment ratio cannot be used 
for the criterion of judgment whether we increase or decrease the temperature. 

As another quantity, we may treat the correlation ratio, the ratio of the corre- 
lation functions with different distances. For an infinite system at the critical 
point, the correlation function decays as a power of r, 

(g(r))^r^ D -^\ (L = oo, f = 0), (4) 
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with the decay exponent rj. Precisely, the distance r is a vector, but we have 
used a simplified notation. Away from the critical point, the ratio of the dis- 
tance r and the correlation length £ plays a role in the scaling of the correlation 
function. Moreover, for finite systems, two length ratios, r/L and L/£, come 
in the scaling form of the correlation function. Then, the ratio of the spin-spin 
correlation functions with different distances r and r' takes the FSS form with 
a single scaling variable, 



if we fix two ratios, r/L and r/r'. In case the correlation length diverges with 
a power law, £ oc t~", Eq. (5) becomes the same form as Eq. (3); however, 
Eq. (5) is also applicable to the case of the KT transition, where the correlation 
length diverges more strongly than the power-law divergence. 

In order to examine the FSS properties of the correlation ratio, we here give 
the result of the 2D 6-state clock model on the square lattice with the peri- 
odic boundary conditions. The 2D g-state clock model is known to exhibit two 
phase transitions of the KT type at 7\ and T 2 (7\ < T 2 ) for q > 4 [20]. We sim- 
ulate the 2D 6-state clock model by the use of the Wang-Landau algorithm [8]. 
We show the temperature dependence of both the moment ratio (m 4 )/(m 2 ) 2 
and the correlation ratio (g(L/2))/(g(L/4)) in Fig. 1. As for the distances r 
and r', we choose L/2 and L/4; we take the horizontal or vertical direction of 
the lattice for the orientation of two sites. For the correlation ratio, the curves 
of different sizes merge in the intermediate KT phase (T\ < T < T 2 ), and 
spray out for the low-temperature ordered and high-temperature disordered 
phases, which is expected from the FSS form of Eq. (5). Then, we can make 
a FSS analysis based on the KT form of the correlation length, 

£ oc exp(c/v / t), (6) 



where t = \T -T KT \/T KT . Using the data of (g(L/2)) / (g(L/4)) for L=12, 16, 
24, 32, 48, and 64, we estimate two KT transition temperatures. The best- 
fitted estimates are Ti=0. 698(4) and T 2 =0. 898(4), in units of J/k B , which are 
compatible with the recent results using the PCC algorithm [12], Ti=0. 7014(11) 
and T 2 =0. 9008(6). On the contrary, as seen from Fig. 1, the corrections to FSS 
are larger for the moment ratio, which makes the FSS analysis difficult. We 
have shown that the correlation ratio is a good estimator especially for the 
KT transition. It is due to the fact that we only use the property of correla- 
tion function; the characteristic of the KT transition is that the correlation 
function shows a power-law decay at all the temperatures of the KT phase. 

We may use the FSS properties of the correlation ratio for the generalization 
of the PCC algorithm. Instead of checking whether the clusters are percolating 
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Fig. 1. The moment ratio (m 4 }/(m 2 ) 2 and the correlation ratio (g(L/2))/(g(L/4)) 
of the 2D 6-state clock model for L = 12, 16, 24, 32, 48, and 64. 

or not, we ask whether the instantaneous correlation ratio g(L/2)/g(L/A) is 
larger or smaller than some fixed value R c . Of course, we can use other sets of 
distances. We decrease (increase) the temperature, if g(L/2)/g(L/4) is smaller 
(larger) than R c . We start the simulation with some temperature. We make 
the amount of the change of temperature, AT, smaller during the simulation; 
in the limit of AT — > 0, the system approaches the canonical ensemble. 

As an example, we apply the generalized scheme of the PCC algorithm to 
the study of the 2D 6-state clock model. We treat the systems with linear 
sizes L = 8, 16, 32, 64, and 128. We start with AT = 0.005, and gradually 
decrease AT to the final value, 0.0001. After 20,000 Monte Carlo sweeps of 
determining T K t(T), we make 10,000 Monte Carlo sweeps to take thermal 
average; we make 20 runs for each size to get better statistics and to evaluate 
the statistical errors. We calculate g(L/2)/g(L/4) to check whether it is larger 
than R c or not. 

We use the FSS analysis based on the KT form of the correlation length, 
Eq. (6). For the L dependence of Tkt(L), we have the relation 

W)=T KT + -^. (7, 

We plot T 2 (L) as a function of l~ 2 with I = InbL for the best- fitted parameters 
in Fig. 2(a). Here, we concentrate on the high-temperature transition T 2 (L). 
The value of R c has been set to be 0.87 and 0.85. The error bars are smaller 
than the size of marks. We should mention that the data with different R c 
are collapsed on a single curve in this plot, which means that b depends on 
R c in Eq. (7) and the difference of R c can be absorbed in the R c dependence 
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Fig. 2. (a) Plot of T 2 (L) and (b) the logarithmic plots of (g(L/2)} at T 2 (L) of the 
2D 6-state clock model for L = 8, 16, 32, 64, and 128, where I = ln&L. The closed 
circles are data for R c = 0.87, and open circles are those for R c = 0.85. 

of b. The present estimate of T 2 is 0.899(3). This value is consistent with the 
estimate of the PCC algorithm using the percolating properties, 0.9008(6) [12]. 

We also estimate the critical exponent i] from the size dependence of the 
correlation at T 2 (L). In Fig. 2(b), we plot {g(L/2)) at T 2 (L) as a function of 
L in a logarithmic scale. For the estimate of r], we use the FSS form including 
small multiplicative logarithmic corrections 

(g(L/2)) = AL-o(\nb'L)- 2r . (8) 

Our estimate is r] = 0.250(3) from the data for R c = 0.87, and r] = 0.265(2) 
from the data for R c = 0.85, which are compatible with the theoretical pre- 
diction, 1/4 (=0.25). 

The detailed description of the generalized scheme of the PCC algorithm based 
on the FSS of the correlation ratio together with its application to the 2D 
quantum XY model of S — 1/2 will be reported elsewhere [21]. 



3 Cluster-flip flat histogram method 

In this section, we consider the combination of the cluster algorithm and the 
extended ensemble method. One calculates the energy DOS, g{E), in the mul- 
ticanonical method [3,4] and the Wang-Landau method [8]; the energy his- 
togram H(E) is checked during the Monte Carlo process. In contrast, the DOS 
for bond number rib, fi(nb), is calculated in the multibondic ensemble method 
[14] or its improvement by Yamaguchi and Kawashima [15]; the histogram for 
bond number, H{rib), is checked in the Monte Carlo process. 

In proposing the broad histogram method, Oliveira et al. [5] paid attention to 
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the number of potential moves, or the number of the possible energy change, 
N(S, E — > E'), for a given state S. The total number of moves is 

J2N(S,E -> E + AE) = N (9) 

AE 

for a single-spin flip process, where N is the number of spins. The energy DOS 
is related to the number of potential moves as 

g(E) (N(S, E -> E')) E = g(E') (N(S f , E> -> E)) E , , (10) 

where (• • ■) E denotes the microcanonical average with fixed E. This relation is 
shown to be valid on general grounds [22,23], and hereafter we call Eq. (10) as 
the broad histogram relation (BHR) for the energy. One may use the number of 
potential moves N(S, E — > E') for the probability of updating states. Alterna- 
tively, one may employ other dynamics which has no relation to N(S, E — > £"), 
but Eq. (10) is used when calculating the energy DOS. 

It is interesting to ask whether there is a relation similar to the BHR, Eq. (10), 
for the bond number. Here, using the cluster (graph) representation, we derive 
the BHR for the bond number. We consider the g-state Potts model as an 
example. With the framework of the dual algorithm [24,25], the partition 
function is expressed in the double summation over state S and graph G as 

Z(T) = J2V (G)A(S,G), (11) 

S,G 

where A(S, G) is a function that takes the value one when S is compatible to 
G and takes the value zero otherwise. A graph consists of a set of bonds. The 
weight for graph G, V (G), is defined as 

V Q (G) = V (n b (G),T) = (e J / fesT - 1)"^ (12) 

for the g-state Potts model, where rib(G) is the number of "active" bonds in G. 
We say a pair is satisfied if Oi = <jj, and unsatisfied otherwise. Satisfied 
pairs become active with a probability p — 1 — e~ J l ksT for given T. 

By taking the summation over S and G with fixing the number of bonds rib, 
the expression for the partition function becomes 

N B 

Z(T)= ]TQK)HK,n (13) 

n b =0 

where N B is the total number of nearest-neighbor pairs in the whole system. 
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Here, is the DOS for the bond number defined as the number of con- 

sistent combinations of graphs and states such that the graph consists of n& 
bonds. Then, the canonical average of a quantity A is calculated by 

{a)t _ E,(^»yo(»,T) (i4) 



where (A) n is the microcanonical average with the fixed bond number rib for 

\ 'ill) 

the quantity A. Thus, if we obtain and (■ • •) during the simulation 

process, we can calculate the canonical average of any quantity. 

In deriving the BHR for the bond number, we follow a method similar to that 
used to derive the BHR for the bond number [22,23]. Instead of using the 
relation between states, we consider the relation between graphs. The number 
of potential moves from the graph with the bond number n b to the graph with 
n b + 1, N(S,G,rib — > rib + 1), for fixed S is equal to that of the number of 
potential moves from the graph with rib + 1 to that with rib, N(S, G', rib + 
1 — > rib). Taking a summation over states S and using the definition of the 
microcanonical average with the fixed bond number rib, we have 

Q(n b ) (N(G, n b -+ n b + 1)L = n(n b + 1) (N(G>, n b + 1 - n b )) n . 1 . (15) 



This is the BHR for the bond number. It should be noted that N(G, rib —> 
rib + 1) is a possible number of bonds to add, and related to the number of 
satisfied pairs for the given state S, n p (S), by N(G, n b — > n b + 1) = n p (S)—n b . 
With use of the microcanonical average with fixed bond number for n v , we 
have the relation 

(N(G,n b ^n b + l)) n = (n p ) n - n b . (16) 



On the other hand, the possible number of bonds to delete, N(G', rib + 1 — > rib), 
is simply given by n b + 1, that is, 

(N(G',n b + l^n b )) nb+1 = n b + l. (17) 



From the BHR for the bond number, Eq. (15), we have 
OK) ."fr 1 (N(G,l^l + l)) nb=l 

n(o) - fi (N(G,l + l^l)) nb=l+1 1 ] 



Then, substituting Eqs. (16) and (17) into Eq. (18), we obtain the bond- 
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number DOS, fl(n b ), as 

/ (n ) — I \ 

\nQ(n b ) = lnfi(O) + £ M^TTi \ ^ 

When calculating the bond-number DOS from the BHR for the bond number, 
we only need the information on (n p ) , the microcanonical average with fixed 
n b of the number of satisfied pairs n p . It is much simpler than the case of the 
BHR formulation for the energy DOS. 

Moreover, in the computation of n p , we can use an improved estimator. If a 
pair of sites belong to the different cluster, this pair is satisfied with a 
probability of 1/Q. If a pair of sites belong to the same cluster, this pair is 
always satisfied. Then, we can employ an improved estimator h p as 

VG)=(l-^)EW) + f. (20) 



where represent a cluster that a site % belongs to. Only the information 

on graph is needed. 

Let us consider the update process for the Monte Carlo simulation. In the 
multibondic ensemble method, a graph is updated by adding or deleting a 
bond for a satisfied pair of sites based on Q(n b ) [14]. We may use the number 
of potential move for the bond number, {N(G, ■ ■ -)) n , for the probability of 
update. Using Eqs. (15), (16), and (17), we get the probability to delete a 
bond, 

(n v ) , + 1 — n b 

P{n b - n b - 1) = * , (21) 

( n p)n b -l + 1 



and the probability to add a bond, 

P(n b n b + 1) = J; + 1 , (22) 

\""P/n b 1 



respectively. 

The actual Monte Carlo procedure is as follows. We start from some state 
(spin configuration) S, and an arbitrary graph G consistent with it. We add 
or delete a bond of satisfied pairs with the probability (21) or (22). After 
making such a process as many as the number of total pairs, N B , we flip every 
cluster with the probability 1/2. And we repeat the process. Since we do not 
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know the exact form of (n p ) n , we use the accumulated average for (n p ) n . The 
dynamics proposed here can be regarded as the flat histogram method for the 
bond number, which we call the cluster-flip flat histogram method. As (n p ) nb 
converges to the exact value, the histogram H{n b ) becomes flat. We calculate 
the bond-number DOS, and then calculate various quantities by Eq. (14). 

As a test, we calculate the L x L x L Ising model on the simple cubic lattice 
with the periodic boundary conditions by using the cluster-flip flat histogram 
method. We show (n p ) n /N B as a function of n b for L = 12 by the solid 
line in Fig. 3(a); we give h^/Nb by the dotted line. The number of Monte 
Carlo sweeps (MCS) is 5 x 10 7 . The difference between the solid and dot- 
ted lines represents the number of potential moves {N{n b — > n b + 1)) /Nb, 
whereas the difference between the dotted line and the horizontal axis rep- 
resents (N(n b — > n b — 1)) /Nb- We note that (n p ) n =0 /N B = 1/2, which is 
expected from Eq. (20). The logarithm of the bond-number DOS, lnQ(n b ), 
obtained by (n p ) nb is shown in Fig. 3(b) as a function of n b . The temperature 
dependence of the specific heat is shown in Fig. 3(c), which reproduces the 
result obtained by the conventional method. 




Fig. 3. (a) (n p ) nb /N B and (b) lnQ(n b ) of the 12 x 12 x 12 Ising model obtained by 
the cluster- flip flat histogram method. The dotted line in (a) denotes n b /NB- (c) 
The temperature dependence of the specific heat per spin, where we have used the 
units of J = ks = 1. 

As another example, we simulate the 3D 3-state Potts model on the simple 
cubic lattice. A first-order phase transition occurs in this model. We show 
(n p ) nb /Nb for L = 12 by the solid line in Fig. 4(a); we give n b by the dot- 
ted line. The number of MCS is 5 x 10 7 . The number of potential moves 
{N(n b — > n b + 1)) /Nb and {N(n b — > n b — 1)) /Nb are given in the same man- 
ner as the case of the Ising model. It is to be noted that (n p ) nb=Q /N B = 1/3 
for the 3-state Potts model. The logarithm of the bond-number DOS, In £l(n b ), 
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obtained by (n p ) n is shown in Fig. 4(b). The temperature dependence of the 
free energy is given in Fig. 4(c). The first-order point where the derivative of 
the free energy has a jump is indicated by the arrow. 




Fig. 4. (a) {n p ) n jN B and (b) lnfi(n 6 ) of the 3D 3-state Potts model of L = 12 
obtained by the cluster-flip flat histogram method. The dotted line in (a) denotes 
Tib/Ns- (c) The temperature dependence of the free energy per spin, where we have 
used the units of J = ks = 1. 

The detailed report on the subject of this section will be published in a separate 
paper [26]. There, the results for the application to the 2D Ising and 10-state 
Potts models will be given, and the efficiency of the method will be discussed. 



4 Summary and discussions 

We have discussed the recent progress in Monte Carlo algorithms. First, we 
have argued the generalization of the PCC algorithm based on the study of the 
FSS property of the correlation ratio, the ratio of the correlation functions with 
different distances. We apply this generalized scheme of the PCC algorithm 
to the 2D 6-state clock model. Since we do not use the percolating property of 
the system, we can apply the PCC algorithm where the mapping to the cluster 
formalism does not exist. It can be applied to many problems. For example, 
the cluster formalism does not work well for frustrated systems, but we can 
use the generalized PCC algorithm. We can also apply the generalized scheme 
to a wide variety of quantum systems. 

Second, we have discussed the combination of the cluster algorithm and the 
extended ensemble method. We have derived the rigorous BHR for the bond 
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number, investigating the cluster (graph) representation of the spin models. 
We have shown that the bond-number DOS fi( n f>) can be calculated in terms 
of (n p ) . We have proposed a Monte Carlo dynamics based on the number of 
potential moves for the bond number, which is regarded as the flat histogram 
method for the bond number. We have applied the cluster-flip flat histogram 
method to the 3D Ising and 3-state Potts models. 
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